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LIST OF SYMBOLS . , 


Symbol 

A = [a.J 
ij 


A* 


I 

MB) 

II A|| 
P R(A) 
P R(A*) 


Definition 

M x N matrix A 

[a ] , an N x M matrix 
J 1 

(Moore-Penrose) generalized inverse of the matrix A 
the identity matrix 

largest eigenvalue of the positive semidefinite matrix B 

X t (AA*) , the spherical matrix norm, or the euclidian vector norm 

AA + , the projection onto the range space of A 

A + A, the projection onto the range space of A* 
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ITERATIVE COMPUTATION OF GENERALIZED INVERSES 
WITH AN APPLICATION TO CMG STEERING LAWS 

’ INTRODUCTION 


The equation 


Ax = b 

will possess a solution if, and only if, b lies in the coordinate space spanned by 
the columns of A, i.e., the range space R(A), and in this case every solution 
can be put in the form 

x g = A + b + (I - A + A) z 

where z is some N x 1 column. That solution with z =0, denoted by x^, has 
the property that for all z, 

IIX. II < II X II 

b - g 

i.e., x^ is the minimum-norm solution. If no solution exists, then there is 

some vector in R(A) which is closest to b in the least-square sense; this vector 
is just x^ = A + b so that for all x 

it b - Ax il < lib - Axil 
b — 

These properties establish the importance of the generalized inverse for 
applications . Unfortunately, there is no truly convenient generally applicable 
formula for its computation. Various elimination algorithms based upon the 
defining. equation, i.e., 

AX A = A 
X A X = X 
(AX)*= AX 
(X A)*= XA 



have been devised, but these are sensitive to the preliminary computation of the 
rank of A. This report describes an iterative scheme for computation of the 
generalized inverse and incorporates the scheme into a FORTRAN subroutine 
which may also be used for computing the inverses of nonsingular square matrices. 


FORTRAN SUBROUTINE DESCRIPTION 


The Fortran subroutine listed below is called by the following statement: 

CALL GNVERS (YO, T, R, M, N, A, BD, AP) 

where A is an M x N matrix; BD is any number not less than ^(AA*) ; and 
AP, the output, is the N x M generalized inverse of A. YO, T, R are matrices 
used by GNVERS and need not be defined by the user. A convenient value for 
BD is 

M 

max y 

1 < i < M A, lb. . J 

1=1 I ill 

where 

b ■ m: = aa * • 

The calling program must contain the following statement: 

DIMENSION Y0( ) , T( ) , R( ) 

The arguments in the dimension statement must not be less than the following 
numbers: N*M, M*M, M*M, respectively. It is assumed that M ^ N. If 
M > N, the generalized inverse of A* should be computed: This stratagem is 
adopted to keep the dimensions of certain intermediate matrices small 
( thereby reducing the number of arithmetical operations and computation 

time) and is acceptable since (A* ) + = (A + )* 

The algorithm implemented via this subroutine forms successive 
iterates Y t , Y 2 , . . . which will converge to A + : The maximum number of 


2 



iterations is denoted by ITERS, and is set to 20 in the listing below. The 
iterations cease when 


max 

E s l < i < N 
1 < j < M 


(n) (n - 1) 

y.. - y.. 

1J 1J 


< TOL (n = 2, 3, 


where TOL is set to 10 in the listing below. If convergence does not occur 
in ITERS iterations, an error statement is printed and the subroutine returns 
the last iterate to the calling program. 


GNVERS makes use of a matrix-multiply subroutine whose calling 
statement is 

CALL FMXMT (M, N, IP, A, B, AB) 

which forms the M x IP product AB of the M x N matrix A and the N x IP 
matrix B. A listing of this subroutine is also given. 


The accuracy of the GNVERS subroutine is illustrated by the following 
example. The 3x4 matrix A was defined as 


A = 


4.604359873-01 6.981586633-01 1.637202877-01 7.932543322-01 
8.176181213-01 -5.241646385-01 9.788190850-01 5.955607116-01 
3.456868410-01 -4.876741769-01 -1.229181702-01 -1.267093084-01 


A careful computation using an electronic desk calculator gave the following 
result: 



7.030320327-01 

5.203275028-01 

-6.275139883-01 

5.241249280-01 


-6.471307646-02 

-2.903730234-01 

8.483851012-01 

1.180262294-01 


1.4932868801-00 
-5.949098340-01 
-1.5520373491-00 
-2. 284458009-02 


J 


while the subroutine GNVERS gave 


+ 

A = 


7.03032032-01 

5.20327502-01 

-6.27513987-01 

5.24124927-01 


-6.47130762-02 

-2.90373023-01 

8.48385101-01 

1.18026230-01 


1.49328688+00 
-5.94909835-01 
-1.55203735 +00 
-2.28445806-02 
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PROGRAM LISTING 


The main subroutine listing is given below. 

SUBROUTINE GNVERS(YO,T,R,M,N,A,BD,AP) 
DIMENSION A( l) , Y0( l) , Yl( l) ,T( l) ,R( l) 

ITERS =20 
TOL= l.E-07 
MN=M*N 
MSQR=M*M 
ALF=1.6/BD 
DOl 1=1, N 
DOl J=1,M 
INDEX1= ( J_l)*N+I 
INDEX2= (l-l)*M+J 
1 Y0(INDEX1)=ALF*A(INDEX2) 

D02 L = l, ITERS 

CALL FMXMT ( M, N, M, A, Y0, R) 

D03 1=1, MSQR 

3 R(I)=-R(I) 

D04 1=1, M 
INDEX1=(J-1)*M+I 

4 R(INDEX1)= 1. (H-R(lNDEXl) 

D05 1=1, MSQR 

5 Y1(I)=R(I) 

D06 1=1, M 
INDEX1=(I-1)*M+I 

6 Yl(lNDEXl)=l. 0+Yl(lNDEXl) 

CALL FMXMT(M,M,M,Y1,R,T) 

D07I=1,M 

INDEX1= (i-l)* M+I 

7 T(INDEX1)=1.0+T(INDEX1) 

CALL FMXMT (N, M, M, Y0, T, Yl) 

E = 0. 0 

D08 1=1, MN 
Y=ABSF(Yl(l) -Y0(l) ) 

IF(Y.GT.E) 9,8 
9 E=Y 

8 CONTINUE 
lF(E.LT.TOL) 99, 10 
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10 D02 1= 1, MN 
2 Y0(l)=Yl(l) 

PRINT 11, ITERS 

11 FORMAT( IX, 28HGNVERS DOES NOT CONVERGE IN, 13, 11H ITERATIONS) 
99 CONTINUE 

RETURN 

END 

The FMXMT listing is 

SUBROUTINE FMXMT ( M, N, IP, A, B, AB) 

DIMENSION A( l) , B( l) , AB( l) 

DOl 1=1, M 
DOl J=1,IP 
L=( J~l)*M+J 
AB(L)= 0. 0 
DOl K=1,N 
Kl=(K-l)*M+I 
K2=( J-l)*N +K 

1 AB(L)=AB(L)+A(Kl)*B(K2) 

RETURN 

END 


DER I VAT I ON OF THE ALGOR ITHM 


The Schulz method for inversion of a nonsingular matrix is discussed 
in Reference 1 and 2. The generalized form of this algorithm, called the 
hyperpower method of degree m, is defined by 


R = I 

- AX 

n 

n 

x 1 = 

X [1 + R 

n + 1 

n \ 


n= 0, 1, 2, 


In Reference 1 it is shown that the hyperpower method of degree 3 is optimum 
in a certain sense, and convergence bounds are given. References 3, 4, and 
5 discuss the hyperpower method of degree 2 (the usual form of Schulz’s 
algorithm) applied to the problem of computing the generalized inverse of a 
matrix A. Convergence bounds and optimal initialization values are investi- 
gated in References 1 and 5. Following the techniques in Reference 4 
Lemma 1 and Corollary 1 establish the hyperpower method of degree m for 
generalized inversion. 
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Lemma 1; The sequence of matrices defined by 


R = P_ / . N - AX 1 

n R(A) n 

x ^ = X /p_ + R + R 2 + ... + R m_1 j 
n+1 n\ R(A) n n n / 

(1) 

n= 0, 1, 2, . . 

(2) 

where 


= A*B q (B^ is some nonsingular M x M matrix) 

, (3) 

X = C A* (C is some nonsingular N x N matrix) 
o o ' o 

(4) 

l|P R(A) - AX o !l< 1 , 

(5) 

and 


l|P R(A*) - V" « 1 

(6) 

converges to A + as n -*■ ® . 

Proof: Since A + is the unique solution of the equations 


“ ’ P R(A) 

(7) 

and 


^ = P R(A*) 

(8) 

it sufficies to prove that 


lim 

'*-' P R(A)-‘ X » ,= 0 

(9) 

lim 

n -“ 1|P R(A*) - V" " ° 

(10) 
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To do this, consider the sequences of matrices 


B = B (v . K . + P + P 2 + . . . + P 
n+1 n \ R( A) n n 


m-l\ 

n / 


c , = C (P_, A . .+ Q + Q 2 + ... + Q m_1 
n+1 n\ R(A) n n n 


where 


p = ~ AA B 
n R(A) n 


Q „ ■ P R( A*) - A * AC 


n 


and (n= 0, 1, 2, . . . ) . Then it is obvious that 
X = A*B 


n n 


( 11 ) 


Moreover, 


X = C A 
n n 


because if we let 


Y = C A* 
n+1 n+1 


and note that 


( 12 ) 


< A * = ®n" 1 ( P R(A*) A * - A * AC n A *) = < C' V ( P R(A) 

‘0 


- AY 


n) 


' - ’ A * /P R(A) - AY n)' 


since 


P„, a * x A* = A*P , A . 
R( A*) R(A) 
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then 



so that Y satisfies equation ( 2) . Since X = Y , X = Y for all n, 
n o o n n 

which establishes equation (12) . 

Since 


AX = P n/1 , AX = AX P . A , 
n R(A) n n R(A) 


we have by an easy induction that 


^ k 

R = P , * 

n R(A) 


k-1 


AX I P , . + R + R + ... + R 
nV R(A) n n n 


(k = 0, 1, 2, . . .) 


m 


Setting k= m gives R n+1 = R r , so that 


P , AX - AX l| < || P / * \ - AX I 
R(A) n+1 R( A) n 


m 


which with equation (5) proves equation (9) . The proof of equation (10) is 
analogous. 

Corollary 1: Let 


Xj ( AA* ) a MAA*) s ... > \ r (AA*) 
denote the nonzero eigenvalues of AA* . If 


0 < “ < xjh?) 

then the sequence defined by 
Y = a A* 


(13) 


(14) 
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R = I - AY 
n n 


( 15 ) 


Y = Y 
n + i 


/i + R + R 2 + . . . + R m A I 
n\ n n n / 


n = 0, 1, 2, ... 


(16) 


converges to A as n — ® . 
Proof: Since 


R(A) 


AA* = AA + AA* = AA* = AA*AA + = AA*P. 


R(A) 


P , . and AA* are commuting hermitian matrices with the same range 
R(A) 

space, so that the eigenvalues of P ^ ^ - AY are 


1 - a X (AA*) 


(i = 1» 2 r) 

0 (i = r + 1, . . . , m) 

hence, by equation (13) 


(17) 


K P B(A) - AY 


o) 


< 1 


(18) 


Similarly 


K P K(A*> - Y o A ) 


< 1 


Now the process ( 2) , begun with equation ( 14) , retains the form 
of equation ( 12) , so 


X / P / . » + R + . . . + R 
n\^ R(A) n 


m -1\ ^ _ m -1\ 

) - X P , v + X(R + ...+R ) 

n / n R(A) n\ n n / 

= C A* AA + + X (r + . . . + R m_1> ) 
n n\ n n / 

= X (i + R + ... + R m_1 ) 
n\. n n / 


Thus, the convergence of equation (16) follows from that of equation (12) . 
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Corollary 2: The process ( 18) is convergent if a satisfies the inequality 


0 < a 


< 


max 
1 < i < 


M 



b. . 
ij 


where B = 


ij 


M 


= AA* 


1 . 


Proof: By Gershgorin’ s theorem 


. t a A *\ mSLX 

^ > s i s i s 


M 


[ 2 ] 

M 

£ 

i=i 


b. . 


Note that it is immaterial whether we take B as AA* or A* A since the 
nonzero eigenvalues of these two matrices are the same. 


APPLICAT ION TO REDUNDANT CMG ASSEMBLIES 


The instantaneous torque output of a single -gimbal CMG is given by 
h = a e x he = hae 

1 1 u 

where a is the instantaneous gimbal rate, h is the CMG momentum magnitude 
(determined by the wheel speed) , and e and e are unit vectors along the 

X u 

gimbal axis and in the direction of the momentum vector, respectively. Then 
e is a unit vector along the instantaneous torque vector h. The torque output 
of N such CMGs is then 



\ + h 2 + 



Vi e i + 


h 2“2 e 2 + 


W 


'N 


where denotes the gimbal rate of the ith CMG and e. now denotes the unit 

vector in the direction of its instantaneous torque vector. This last equation 
can be put in the following form: 

h = Co- 
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where C is the 3 x N matrix whose columns are h^e^, h^e^, * * * ’ ^N 6 N an< ^ 
a is the column of gimbal rates. Then, we have that 1 ' 1 


B = A* A = 

i 2 ji. 

h i e i e i 

h l h 2 e i* 6 2 • • • 

h lVl* e N 


h h e * e 
2 12 1 

' • • 

h 2 h N e 2* C N 


VA* e i 

V2 e N* e 2 ' ' • 

\,V e N _ 


MA*A) 


N 


max , , 

1 < i< N H i j 

3=1 


e. ©. 

i J 


max 

1 < i < N 


N 

Y h. h. :£ Nh 2 

& 1 3 


since | e.*e. | 
i J 

used. 


1 for each i and j; h is the nominal momentum of the CMGs 


If T is a commanded torque output, we wish to ensure that h = T by 
c c 

choosing a; the procedure by which this choice is made is called a steering 

law. One such law, which is optimum in a certain sense, is 

a = C + T 

c 

While this steering law has exhibited satisfactory performance in a 
number of simulations, a number of complicating factors must be considered 
before the law can be considered for implementation. Discussion of these, 
which have mainly to do with situations in which the columns e^ become 

coplanar (a hangup condition), is beyond the scope of this report. 
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